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1. INTRODUCTION 


The GEOS-3 investigations at UCLA originally proposed a 
mixture of geoid determination and geophysical analysis and 
interpretation. While efforts in all these directions were 
undertaken, it has become apparent that, taking into account 
work elsewhere, the area in which we could most significantly 
contribute is the analysis and interpretation, which is there- 
fore emphasized in this report. 

Variations in the external gravity field — expressed either 
as geoid undulations or gravity anomalies — indicate, of course, 
lateral variations in density within the solid earth. In any 
laterally coherent material, such density variations must entail 
shear stresses. The finite strength of rocks limits these 
stresses; there will always be fracture and flow until the 
density variations are removed or compensated in some way. 

Because (1) this flow has a rate highly non-linear in stress; 

(2) the rheological properties of rocks are 3trongly tempera- 
ture dependent; and (3) the stress for a given amplitude of 
density irregularity is proportionate to the wavelength, we 
should expect even the relatively plastic oceanic lithosphere 
and upper mantle to show a gravity signal whose broad wavelengths 
correspond to isostatic compensation of the rise to basin topo- 
graphic drop around a3thenospheric depth, but whose local varia- 
tions indicate departures from this overall "t* 5 " sweep. This 
oceanic isostasy differs from the more familiar continental 
variety in that iu is mainly thermal, rather than compositional, 
in origin: in general, there is a positive correlation of temper- 

ature with topography, hence a negative correlation of density 
with topography, such that the total mass of a column of matter 
down to somewhere around the asthenosphere is about the same 
everywhere. Models of thermal convection define the distribu- 
tion of this isostatic compensation. Since, in principle, 
pressure effects can result in differences from this equal 
loading, perhaps the term 'isodynamics' should be used rather 
than 'isostasy'. But these differences are small in models 
developed so far. 

Given the foregoing considerations, together with the 
likelihood that features deeper in the mantle not connected 
ir an obvious way to plate tectonics are contributing to the 
gravity signal, any rational procedure of analyzing or inter- 
preting the ocean geoid and topography will take into account 
their spectra. Procedures which have been applied can be 
classified are varying mixtures of four archetypes: (1) 

statistical analyses, in which auto- and cross-covariances and 
their spectral transforms, or admittances, are determined from 
gravimetry and bathymetry (e.g., McKenzie & Bowin, 1976; Kaula, 
1977); (2) removal of empirically determined mean trends and 

mapping of residuals (e.g., Sclater et al., 1975; Cochran & 
Talwani, 1977); (3. iirect modelling ir which the effects of 
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a specified thermal and mechanical model on observables are 
calculated (e.g., McKenzie et a]., 1974; Schubert et al., 1976); 
and (4) inverse modelling, in which the best values for the 
parameters of a specified model are inferred from data (e.g., 
Burkhard & Jackson, 1977). At present levels of data coverage, 
technique, and insight, any particular study is a trade-off 
between the different aspects: elaborate statistical analyses 

assume rather simple interior models, while complex physical 
models account for only certain features of the data. We 
discuss here efforts in categories (1) and (3) in turn: 
statistical analyses and 1 ithospher ic -asthenospher ic modelling. 
As will be seen, the necessity of removing effects of orbital 
errors from altimetry make statistical analyses willy-nilly 
incorporate (2), removal of empirically determined trends. 


2. STATISTICAL ANALYSES 
2 . 1 Admittance Analysis Fundamentals 

McKenzie (1967) suggested that the relationship between 
gravity (or geoid) anomalies and topography should be analyzed 
in terms of the admittance. If the gravity anomalies Ag and 
ocean depths b along a track are expressed in terms of complex 
Fourier series, 

Ag * Jg (k) e i,CX "] 
b- | • 

where x is along-track location and k is the wave number, 
then, if there are no other density irregularities than the 
topography, the complex admittance Z(.1is defined by (McKenzie 
& Bowin, 1976) : 

g ( k ) = Z ( k ) b ( k ) 1 (2.2) 

The first ta.sk is to estimate the admittance from the gravity 
and bathymetry over specific tracks, and the second to in- 
vestigate what models produce admittances most closely approxi- 
mating the observed admittance. 

To estimate the admittance from N equal blocks of data, 
the mean and trend from each block r are removed by least 
squares, the Fourier transforms g r (k), b r (k) obtained, and the 
auto- and cross-power spectra determined from 


s k (gg 0 

, N 

g r (k)g*(k) ~j 

S k (bb) 

S R (gb) 
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where the asterik denotes the complex conjugate. The ad- 
mittance is then obtained from 

Z (k) - S k (gb)/S k (bb) (2.4) 

2 

The coherence y (k) is defined by 

y 2 (k) - S k (gb)S* (gb)/S k (gg)S k (bb) (2.5) 

All these functions can equally well be determined for geoid 
heights N as for gravity anomalies Ag, of course. 

2 . 2 Gsoid and Seafloor Ba t hymetry 

We utilized the sea surface heights above the ellipsoid, 
corrected for tidal variations, produced by NASA Wallops Flight 
Center. A standard deviation of 0.1-C.3 m was commonly found 
for individual altitude measurements, although values of over 
1 m may be found on bad passes. Measurements tended to be 
poorer ovei very calm seas or upon approaching land. This 
standard deviation of 0.1-0. 3 m is probably representative of 
the relative accuracy of the measurements along an individual 
satellite pass. Orbit errors were indicated by large residuals 
at track crossings. From our work, the orbital errors appeared 
to have wavelengths of 10,000 km or more. The estimated re- 
lative standard deviations of 0.1-0. 3 m is certainly small 
enough to resolve 3-6 m perturbations in the geoid along long 
tracks. Some problems may be experienced in determining areal 
geoids involving numerous tracks, however, if these long wave- 
length orbit analysis errors are not reduced. 

The locations of the satellite tracks used in this study 
are shown in Figure 1 along with their GEOS-3 orbit number. 

The data was selected from the Pacific Ocean to avoid the 
coiplications of isolated volcanic centers. (Analyses of such 
areas will be undertaken when more complete data coverage is 
available). Because of the limited current coverage, emphasis 
thus far has been placed upon analysing tracks oriented as 
closely as possible to perpendicular to the direction of 
seafloor spreading. The available tracks entailed considerable 
compromise from ideal, however. 

The geoid undulations are presented above the reduced 
seafloor bathymetry in Figure 2. Bathymetry for the North 
Pacific is from Chase et al. (1970), and for the South Pacific 
from Mammerickx et al. (1975). Depths have been corrected 
for sound velocity (Matthews, 1939) and for sediment loading 
(Ewing et al., 1968; Ewing et al., 1969). Isochron locations 
used to calculate residual depths (not shown) are taken from 
Heron (1972), Larson and Chase (1972), Larson and Hilde (1975), 
and Hilde et al. (1976), and the empirical depth versus age 
relation is from Parsons and Sclater (1977) . 
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Figure 1. Locations of altimetry tracks. Ocean rises: ++♦+++; 

subduct iori zones: iiinnni . 
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Before any admittance analysis is undertaken, it is worth 
examining the qualitative correlation of the geoid with either 
the bathymetry or the residual depth. The profiles in Figure 2 
reemphasize the results of Cochran and Taiwan! (1977) s the 
intermediate wavelength geoid (200-1000 km) shows reasonable 
correlation with most of the volcanic features, but the longer 
wavelength geoid (greater than 1000-2000 km) shows only slight 
correlation with the topography (Figure 2) . No improvement 
is realized in the visual correlation at long wavelengths by 
additional smoothing of the profiles. 

Strong positive correlations between geoid and topography 
can be seen where track 1993 crosses several seamounts near 
the Chinook Fracture Zone and where it intersects the Emperor 
Seamounts. A positive correlation is also noted where track 
1962 crosses the Line Island Ridge (at about 4500 km) ; almost 
no geoid anomaly exists where the track crosses the Cook 
Islands implying more local compensation. The tracks over the 
youngest lithosphere (737 and 752) exhibit only low amplitude 
and long wavelength geoid anomalies as do tracks 1243 and 319 
also over young lithosphere. Except for track 752 which showns 
a strong negative correlation between topography and geoid, 
the topography over this young lithosphere appears to be 
compensated at shallow depths. Some evidence for a negative 
correlation at long wavelengths is also observed in tracks 
1993 and 1962 away from the volcanic features. The evidence 
is not definitive, however. 

2 . 3 Application of Admittance Analysis 

2.3.1 Isolated Track Analysis 

The geoid height and bathymetry profiles were Fourier 
transformed, and admittances (eq. 4) and coherences (eq. 5) 
calculated for each satellite track using a band-averaged cross- 
spectral technique. The results are shown in Figure 3, with 
the complex component translated into phase. Tor the band- 
averaging used (Q * 0.25) and 95% significance level the co- 
herence is about 0.6. As can easily be seen, few of the 
coherences between geoid and topography exceed this value 
indicating that less than 33% of the energy is coherent between 
the two profiles. The overall lack of coherence between the 
geoid and topography is further accentuated by the large scatter 
in phases. Simple causative relations between the two would 
result in either a 0° or 180° phase. Calculation of the ad- 
mittance between geoid and residual depth produces comparable 
or lower coherences and equally poor phases. The amplitude of 
the admittance is indistinguishable from those calculated from 
the corrected bathymetry. This is largely because ill the 
tracks have been detrended before processing using a second 
order polynomial. The second order polynomial removes most 
of the longest wavelength variations associated with the age 
dependence. It should also be noted that most of the tracks 
are nearly parallel to the isochrons. It is interesting to note 
that the coherences that are equal to or greater than 0.6 are 
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Figure 3. Results of the spectral analysis. The 95« significance level 
for the coherences is approximately 0.6. The error bars on 
the amplitudes represent one standard deviation. No error 
bars are shown for the phases because at these low coherences 
one standard deviation is nearly +180°. Not that there is 
usually a large change in phase between the wavelengths of 500 
and 1000 km. (-*) 1993; (x) 1962; (A) 319; (D) 1243; (o) 737; 
(V) 752. 
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usually for wavelengths greater than 1000 km (see tracks 1962 
and 752). The high coherence for track 1962 may result from 
the dominance of the Line Island Ridge in the profile; no 
volcanic features stand out on track 752. The significance 
of this result is again somewhat diminished by the scatter in 
the phases. 

In an effort to obtain more statistically significant 
admittances, the results from the three oldest tracks, 1993, 

1962 and 2488 have been combined. The results are shown in 
Figure 4. The averaging has apparently caused little change 
in the admittances. The coherences are still statistically 
insignificant for wavelengths smaller than 1000 km; however, the 
coherences for wavelengths greater than 1000 km have now become 
more noteworthy because the 95% significance level has dropped 
to 0.42. If the low coherences found were simply the result of 
a high noise level, increasing the number of samples by averaging 
should improve the coherence. The lack of improvement may imply 
that the physical mechanism /oiponsible for the geoidal undula- 
tions of 300-1000 km wavelength are truly unrelated, or &t least 
only weakly related, to seafloor topography. The non-zero phases 
also appear to be unchanged by the averaging and do not manifest 
the instability often associated with high noise levels. A 
phase reversal may still be surmised from the rapid change of 
phase near 1000 km wavelength. 

Of the three parameters calculated, the phase is by far 
the most uncertain. Error bars at the 95% confidence level 
span the entire 360° cycle and are thus not shown. The phase 
is most sensitive to the details of the bathymetry and the 
sediment corrections. In areas such as the South Pacific, it 
is not known how accurate the bathymetry and isopachs are, even 
at long wavelengths. Poor phase estimates are also found in 
the NE Pacific where both are well known suggesting that the 
origin of the phase uncertainties should be sought elsewhere. 

One advantage of shipboard gravity measurements is that the 
bathymetry can be concurrently measured; however, the short 
wavelength features dominate the signal. With these points 
in mind, it is difficult to place even a positive or negative 
sign on the admittance at long wavelengths. However, the 
trend in most of the phase curves suggests that theie may be 
a phase reversal somewhere between the wavelengths of 500-1000 km. 

2.3.2 Discussion 


The low coherences and small admittances at intermediate 
wavelengths (200-1000 km) support the conclusion that the 
intermediate wavelength topography may be compensated at shallow 
depth (Figure 4). In fact, both this study's results and those 
of McKenzie and Bowin (1976) might be explained by compensation 
in a layer of varying thickness (Airy compensation) . McKenzie 



and Bowin rejected this model because the average thickness of 
this layer must be much thicker than the oceanic crust. If 
this Airy compensation model accounted for all the observations, 
a much higher coherence would be expected and the admittance 
should be real. 

The model favored by McKenzie and Bowin was their plate 
model. In this model the intermediate and long wavelength 
topography is compensated by bending an elastic plate of total 
thickness 2h (solid lines in Figure 4); the short wavelength 
topography is not compensated. The plate model reproduces 
the observations with wavelengths less than 1000 km rather well 
if the plate half-thickness (h) is taken to be 5-10 km and the 
shear modulus between 10^° and 10^1 Nm“2. Again, however, the 
consistently low coherences argue against an admittance wholly 
determined by the plate model. 

The above models require a high coherence between geoid and 
topography because the source of both lies within the litho- 
sphere. The low coherences might indicate that the source of 
the geoid undulations lie below the lithosphere. The gravity 
anomalies and deformations arising from convective processes 
result in a real or zero admittance (McXenzie, 1977). The zero 
admittance is the result if large viscosity variations completely 
decouple the convective stresses from the overlying plate. The 
decoupling effectively reduces or eliminates the coherent surface 
deformation. This appears to be the case best represented by 
the data. Not enough data currently exists to put rigorous 
bounds on the nature of the decoupling zone, however. 

McKenzie (1977; Appendix C) gives admittances for convection 
partially decoupled from the lithosphere by a low-viscosity zone 
of thickness d. The admittance depends upon the thickness of 
the layer (d) and the viscosity contrast between the convecting 
region and the low-viscosity layer. The thickness of the low- 
viscosity layer is determined from the wavelength at which the 
admittance changes phase (1000 km for the averaged data in 
Figure 4) and the viscosity contrast by both the range of wave- 
lengths over which the admittance remains negative and by the 
largest negative amplitude. Unfortunately, we have little 
control over the viscosity contrast due to the lack of longer 
wavelength data, but the thickness of the low-viscosity layer 
is estimated to be 200-300 km. There is no resolvable age 
dependence, although the profiles over younger lithosphere 
change phase at shorter wavelengths. 

Kaula's (1977) statistical analysis of gravity and topo- 
graphy over oceanic regions agrees well with the results pre- 
sented here. Similar negative admittances are implied for 
wavelengths shorter than 5000 km. (Kaula's analysis only 
included harmonic degrees as high as 36, so no lower limit can 
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Figure 4. Several models are superimposed on the averaged admit- 
tances. Short wavelength admittances from McKenzie and 
Bowin (1976) are shown as 3olid rectangles. The solid 
curves are from the plate model; the upper curve is for 
h = 10 !cm and the lower is for h * 5 km. All other 
parameters are as in McKenzie and Bowin. The dashed 
curve is the decoupling model of McKenzie (1977) for 
M = 100 and d = 200 km. 
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be estimated.) However, combining this 5000 km upper limit 
with the lower limit of 1000 km from this study suggests a 
viscosity contrast of about 100 between the asthenosphf re 
and the mantle below. 

A decoupling model also admits an explanation for the non- 
real admittances. Even though the convective stresses create 
negligible seafloor topography, topography can result from 
thermal expansion. Since the lithospheric plate may move re- 
lative to the convecting mantle, the finite thermal conductivity 
of the lithosphere will cause any topography created by thermal 
expansion to be out of phase with respect to the mantle heat 
sources responsible for the geoidal undulations. If this is 
the case, some estimate of the relative motion between the 
lithosphere and convecting mantle might be obtained by studying 
the phase differences between geoid and gravity and topography. 

It is not entirely clear whether a single physical process 
correctly describes the data. For instance, the longer wave- 
lengths (greater than 1000 km) may be controlled by the con- 
vection and the viscous decoupling, but the intermediate wave- 
lengths may result from compensation of inhomogeneities within 
the plate through a plate model (McKenzie and Bowin, 1976) or 
Airy model. The intermediate wavelength admittances do fit a 
plate model within reason. Or, both processes may simultaneously 
occur. 

Further speculations on the origin of the geoid and sea- 
floor topography r-'ist await more comprehensive data. Thus far 
only two-dimensioi. 1 admittance have been studied and not neces- 
sarily in the optimum orientations nor o.er a suitable spread of 
wavelengths. Consideration of the geoid over a large oceanic 
area will belter enable us to relate the long wavelength geoid 
to the regional tectonics and to appreciate t^ a contribution 
ot "isolated" volcanic features to the regional geoid. Admit- 
tances from three-dimensional convective systems it, ay result in 
totally new interpretations. 

2.3.3 Areal Admittances over the Hawaiian Swell 

To further investigate the relationship between geoid 
and seafloor topography two-dimensional admittances were 
calculated for the area shown in Figure 5. Admittances band- 
averaged according to wave number were nearly identical to 
those found in the track analysis. The coherences were also 
found to be similarly low. 




Figure 5. Study area. Grey dep 
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Directional band-averaging, however, resulted in anisotropic 
admittances and coherences. Admittances crossing east-west 
striking structures (from north-south profiles) were smaller 
but more coherent than those from east-west profiles. Since 
the geological structures in this area predominantly strike 
east-west, this is further evidence for a coherent geoid only 
ov .r volcanic structures. Unfortunately, the area considered 
for analysis was too small to provide good estimates of the 
long wavelength admittances. 

2.3.4 G lobal Statistical Analyses 

Covariance analyses can also be carried out on a sphere: 

C (gb, T ) = --^// 2n g(r)b(T)dAdr, (2.6) 

8Il*To 


where dr is an element of area on a unit sphere of total surface, 

T, A is azimuth, and t is arc distance, and Legendre polynomial 
transforms (analogous to Fourier transforms) applied to obtain a power 
spectrum: 


igb) 


/3T+T 

T~ 


P io (cos 


x)C(gb, T)sintdT 


(2.7) 


where l is spherical harmonic degree and P. is the zonal 
harmonic normal i zed to mean square unity. Tfie estimated ad- 
mittance is then 

Z l = °J(gb)/o 2 (bb) . (2.8) 

Such analyses were applied to global data sets of 5°>:5 0 
means of gravity anomalies based on the NASA-Goddard SFC 
GEM-6 solution (Smith et al., 1976) plus DoD l°xl° means, and 
of topography compiled by Lee (Balmino et al . , 1973). These 
data sets were divided into 15 terrain types. Also included 
in he analysis were plate velocities and boundary types. 

The details are described in Kaula (1977). Given here are 
the results for the combination of oceanic terrain types. 
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6. 
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< 16-18> 

2. 
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0.01 
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1.3 

0.004 

-0.01 

-0.002 

<25-36> 

1.7 
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0.00 

0.000 

Sums 

219. 
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The low admittances in the 7-36 degree range — corresponding 
to wavelengths of 6000 down to 1100 km — are consistent with the 
findings from the altimetry profiles described above. The larger 
positive admittances at 3 to 6 are evidently the consequence of 
seafloor spreading, which trend has been removed from the track 
data discussed above or by McKenzie & Bowin (1976). The s: an 
reversal between degrees 2 and 3 mimics what McKenzie (197V, 

Fig. 6a) predicts for a 3000 km deep convecting layer of non- 
deformable lower boundary, but this seems too vertiginous a 
leap between model and data to take seriously. 

It is planned to redo the global covariance analyses 
utilizing some improvements in technique (Kaula, 1978). 


3. LITHOSPHERIC THINNING 

In this study we will analyze a region in the Central 
Pacific containing a large portion of the Hawaiian chain. Our 
aim will be to show that the major portion of the Hawaiian swell 
is caused by a local thinning of a higher density lithosphere 
on top of a lower density asthenosphere . To accomplish this, 
the relationship between geoid heights and seafloor topography 
will quantitatively be analyzed at intermediate wavelengths 
(300<X<3000 km). Cochran and Talwani (1977) used gravity data 
to qualitatively determine this relationship over many oceanic 
areas. They found that there is generally a poor correlation 
except over active spreading zones or mid-plate volcanic areas. 

A positive correlation between "tred air’gravity and topography 



Figure 6. Seafloor topography averaged into areas 203 km on a side 
with best fitting plane removed. Units ate decameters. 
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or geoid height and topography can be produced by three 
mechanisms: support of topography by stresses within the 

lithosphere, dynamic support at the base of the lithosphere, 
undulating lithosphere/asthenosphere boundary. McKenzie (1967) 
has shown that the first possibility can be discounted for 
wavelengths greater than 1000 km. Cochran and Talwani (1977) 
have attributed the correlation associated with spreading 
ridges to cooling and subsidence of the lithosphere as it 
moves away from the ridge crest. We feel that this thermal 
mechanism also produces the observed correlation over the 
Hawaiian chain. 

In the first part of the analysis we will determine that 
part of the ocean geoid which is coherent and in phase with 
the topograpl . A local compensation mechanism will then be 
proposed to reproduce this part of the geoid. Out of the many 
possible local compensation mechanisms only Aivy-Heiskanen 
compensation using two compensation depths will be investigated. 
Short wavelength topography will be compensated within the upper 
part of the lithosphere and longer wavelengths will be compen- 
sated at the lithosphere/asthenosphere boundary. This lodel 
ir lies a variable thickness lithosphere. Since geoid heights 
alone cannot distinguish between different compensation mechanisms, 
we will try to justify our choice by showing that the variable 
thickness lithosphere model can explain the subsidence of the 
islands along the Hawaiian-Emperor chain. 

Airy compensation for the Hawaiian swell requires that the 
lithosphere beneath the swell be tens of kilometers thinner than 
the surrounding lithosphere. As pointed out by Detrick and 
Crough (1978) this thinning must occur over a time interval of 
less than 5 my. The characteristic time for a heat perturbation 
to conduct through 30 km of lithosphere is about 30 my. Unless 
very large heat anomalies are introduced, the lithosphere cannot 
be thinned fast enough by a purely conductive mechanism (Crough 
and Thompson, 1976) . A way to produce rapid thinning is to 
advect heat to the base of the lithosphere and let the lithosphere/ 
asthenosphere boundary migrate in time. We propose a simple 
one-dimensional heat flow model with a moving solid/liquid 
phase boundary to account for both the rapid thinning and later 
gradual thickening of the lithosphere along the Hawaiian chain. 

3.1 Data 


The region of interest is shown in rigure 5. The area 
outlined is a 3248 km square enclosing much of the Hawaiian 
swell. For analysis purposes the area is divided into 256 
•c^-ual equates- of 2(13 km on $ side (approximately 2°x2°) . Geoxd 
and bathymetric averages are estimated for each square. A 
mean and best fitting plane are removed from both before a two 
square cosine taper is added to the perimeter of the region. 
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Bathymetry . Seafloor bathymetry was taken from Chase 
et al~ (1970) and was corrected for water velocity using 
Matthews tables (Matthews, 1939). No sediment corrections 
were applied to the data because the sedimentary cover is 
regionally uniform (Ewing et al. , 1968) and the large averaging 
areas largely negate the influences of local accumulations. 

Figure 6 shows the results after averaging. The WNW striking 
island chain can clearly be seen. The deep and rise flanking 
the island chain do not show up after averaging because the 
box size is greater than the wavelength of these features. 

The Hawaiian swell shows up clearly with a strike similar to 
the island chain and a half-width of about 1000 km. The swell 
is elongated on the north but drops sharply to the south at the 
southeast end of the chain. The mean depth of the region (5169 
m) fits the empirical depth-age curve at about 60 my; the 
magnetic age of the area is closer to 100 my. 

Geoid . GEOS- 3 altimeter data was averaged into the same 
squares as the bathymetry (Figure 7) . The dominant feature 
in the geoid is the long wavelength high over the Hawaiian 
swell region with an amplitude of roughly 10 m. The central 
portion of this high has a feature elongated in the WNW direction 
This feature probably reflects the part of the geoid produced by 
the island chain. The remainder of the high is roughly circular 
with a half-width of about 1200 km in the NS direction and 
about 1600 km in the EW direction. The peak- to- trough amplitude 
in the EW direction is about 12 m and in the NS direction is 
only 10 m. 

Geoid data was used in this analysis for a number of reasons 
First this data is best suited for studies at intermediate wave- 
lengths. At the shorter wavelengths undulations contain very little newer. 

AL the longer wavelengths (angular order <10) the GEOS-3 derived 
geoid may contain inaccuracies due to satellite tracking errors. 
The second reason this data was used is because Watts (1976) 
has used surface gravity data to infer the lithospheric and 
upper mantle density structure over a region roughly the same 
as that shown in Figure 5. Since the GEOS-3 geoid and the 
surface gravity measurements are independent data sets, it is 
interesting to compare the results of the two analyses. 

3. 2 Airy Compensation Models 

There are a number of isostatic compensation models one 
can use to calculate the geoid from the topography. Walcott 
(1969) used a flexure model to fit the observed gravity anomalies 
over the Hawaiian region up to 300 km away from the islands. 

The best fitting flexural parameter was found to be 110 km. 

This model of compensation is important at wavelengths near 
100 km but is less important at longer wavelengths. 

Other important compensation mechanisms include Pratt and 
Airy-Heiskanen compensation. For these types of compensation 
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Figure 7. Ceold heights averaged Into areas 203 kin on a side with 
best fitting plane removed (meters). 
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the mass in every vertical column is the same from the surface 
to the depth of compensation. For the Pratt model this is ac- 
complished by varying the density in adjacent columns. The 
Airy model uses u constant density layer of variable thickness. 
Local compensation can be composed of any number of layers and 
any combination of the Pratt and Airy models. 

The geoid for both of these models can easily be calculated 
using a method proposed by Parker (1972). Let z be positive 
upward and place the origin at a depth z„ below the sea surface. 
z Q is chosen so that the maximum undulation in the density 
interface (i.e., seafloor topography) is minimized. As an 
example, take the calculation of the geoid height N from un- 
compensated seafloor topography h(x,y). Let F(1 denote the 
Fourier transform of the quantity in brackets. Let k, p and 
Pw be the wavevector, cructal density and sea water density 
respectively. Then, 

2ttG(p -p ) v 00 I £1 n-2 

F[N] * — --exp (- | k | z Q ) l -LJ- F[h n (x,y) ) (3.1) 

where g - mean acceleration of gravity. 

These are the exact equations for N due to the seafloor- 
topography. The series converge as long as (max ( h 2 (x,y) ) ) 1 

Fcr our case this means the seafloor topography can not be 
above the sea surface. The smaller the above ratio is made 
the faster the series converges. For all of the models cal- 
culated this ratio was less than 0.5. Convergence of the series 
was more than adequate after summing six terms in the series. 

The geoid height for an Airy model was calculated by super- 
imposing several surfaces of density contrast at different 
depths. Suppose Nj_(x,y) is the geoid anomaly due to the un- 
compensated topography. Let p^ be the density of the crust 
and P 2 the density of the lithosphere. If z^ is the mean 
depth of compensation, then the Fourier transform of the geoid 
height is given by 


FIN] = FINjl+^lPj-Pjlexpl-ISIz^ 


p. -p. 


— ^h (x,y) } n ] 
"°1 

(3.2) 


This can be generalized to any number of layers with different 
wavelengths of topography being compensated at different depths. 

* 9 • • 

The model adopted to satisfy the geoid data is shown in 
Figure 8. Short wavelength topography (less than 800 km wave- 
length) is compensated at crustal depths while longer wavelength 



Diagram of Airy lsostatic model. h(x,y), hjtx.v), and r^lx.v) are measured positive 
upward from zq. zj and z, respectively. lt(x,y) is th.' seafloor topography. hj(x,v) 
and h2(x,y) are calculated from h(x,v) is the wavelength domain. 

y.*) “ H( ^l ~ 800 )h( '' ) (p -P ) • h 2 ( ^ ) " H( 800 " -p f ) 

* • • • ' C ^ ^ ^ ^ 1 

H( ) Is the heavyside step function. * 9* • m 
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features are compensated at the lithosphere/asthenosphere 
boundary. The actual depth of compensation is most likely 
a much more complicated function of wavelength. A model 
with one depth of compensation is the simplest. This model 
was tested for a variety of compensation depths. These models 
fit poorly to the observed geoid at either the long or short 
wavelengths depending upon the depth of compensation. Because 
of this it was necessary to use a model with two depths of 
compensation. The crustal density (p c ) and mean thickness (z^) 
were taken from Hussong (1972) and were not modelled. The only 
free parameters in the model were the mean depth of compensation 
for the long wavelength topography (zo) and the lithosphere/ 
asthenosphere density contrast (p ff -p ). 

* <31 

Before discussing the results of the modelling it is worth- 
while commenting upon what exactly it is that the models are 
to match. An Airy model requires that the geoid be nearly 
coherent and in phase with the topography. It is therefore 
unreasonable to expect any model to account for all of the 
observed geoid. In order to determine what portion of the 
geoid can be modelled, the geoid was separated into two parts: 
one coherent and in phase with the topography, and one inco- 
herent or out-of-phase with the topography. This separation 
was accomplished by calculating the band-averaged two-dimensional 
transfer functions between the topography and geoid. The coherent 
and in phase part of the geoid is reconstructed from the topo- 
graphy and the real part of the transfer functions. It should 
be noticed that the part of the geoid that is incoherent or 
out-of-phase with the topography (Figure 9) accounts for about 
60% of the total energy. The rms residual from no model can 
be smaller than the rms standard deviation (2.71 m) of this 
part of the geoid. 

A comparison of figures 7 and 9 reveals that there is a 
large portion of the Hawaiian geoid anomaly that is not reflected 
in the topography. It is important that this unrelated part 
of the geoid not be included in any regression of geoid and 
topography. Regression slopes obtained without removing this 
incoherent geoid will be strongly biased and will indicate 
deeper compensation depths. 

The best fitting models found have mean depths of compensa- 
tion between 50 and 90 km. This range and the best fitting 
model were determined by two methods. First we subtracted the 
model from the observed geoid and calculated a rms deviation. 
Taking the depth of compensation at 10,000 km (uncompensated 
model) we found an rms deviation of 12.9 m. Figure 10 shows 
this model. It looks very similar to the averaged seafloor 
topography (Figure 6) . By comparing this model to the observed 
geoid (Figure 7) one can see that the amplitudes are much too 
large". The best f it^ing • comperes-a ted model foond by this ^c>hod 
has a Z 2 of 70 km and an rms deviation of 3.04 m. Figure 11 
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Figure 10. Geold height calculated fr^m uncompensated seafloor topography 
with density of ?.87 gm/cnr. 1’nlts are In meters. 
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shows the observed geoid minus the model for zj ■ 70 km. 

Using this method we could not determine a range for the depth 
of compensation because the incoherent and out-of-phase part 
of the observed geoid produced a major part of the rms deviation. 
By fitting models to the in phase and coherent part of the 
geoid limits could be part on Z 2 > Again, a best fitting model 
with Z 2 * 70 km was found with an rms deviation of 1.2"* m. 

Depths of compensation of 50 and 90 km have rms residuals of 
1.45 and 1.48 u respectively. The Table shows the rms deviation 
between, the model and observed geoid, and model and the coherent 
and in phase geoid, for a number of compensation depths. 


Depth (km) 

Observed Geoid 
-Airy Model, Meters 

Coherent 

-Airy 

and in Phase Geoid 
Model, Meters 

50 

3. 12 


1.45 

70 

3.04 


i. 27 

90 

3. 10 


1. 48 

110 

3.17 


2.43 

lj 0 

3.49 



10 5 

12.97 




The density contrast (Pj-pJ is the parameter that determines 
the amplitude of the undulations of the lithosphere/asthenosphere 
boundary needed to support the surface topography. Models with 
density contrasts of 0.06 to 0.12 gm/cm 3 produced nearly identical 
geoids. Because of this insensitivity, it was not possible to 
place definitive bounds on the density contrast. A lower bound 
might be estimated by considering the maximum thinning possible. 
The i.iry model predicts that the lithosphere/asthenosphere 
boundary will intersect the crust/lithosphere boundary if the 
density contrast is less than 0.06 gm/cm 3 . A density contrast 
'f 0.1 gm/cm 3 gives a more reasonable maximum thinning of 35 km. 
This density contrast produces a Lithospheric thickness beneath 
the swell corresponding to a thermal age of about 35 my (Parsons 
and Sclater, 1977) . Figure 12 shows the two compensation 
boundaries for tracks along and perpendicular to the island 
chain for a density contrast of 0.1 gm/cm 3 . The location of 
these tracks are indicated on Figure 6. The thickening is ob- 
served after the volcanically active islands are passed. The 
maximum thinning does not appear to take place directly under 
Hawaii but slightly to the west of it. 

3 . j T hermal Model 

The empirical depth/age and heat/flow age relations can 
be fit by a thermal spreading model with a moving phase boundary 
and heat flow added from below (Oldenburg, 1975). This same 
model can be used to explain the rapid uplift of the seafloor 
topography at the southeast end of the Hawaiian chain and also 
the gradual subsidence of the older islands. 



DEPTH (km) 


26 



Figure 12. Cross sections of the Hawaiian chain assuming Alrv compensation (Fig. 8). Tracks 
A-B and C-D are shown on Fig. 6. In each cros_ section, the top dotted line Is 
the mean ocean depth, the top solid line Is the crust /l lthosphere interface, and 
the bottom solid line Is the 1 lthosphere/asthenosphere Interface. 

Within the solid lithosphere the temperature is governed 
by the usual heat conduction equation 


3T 

Tt 



(3.3) 


with an upper boundary condition of T(0,t) = 0 and T(Z(t),t) = 
at the lower boundary. The bottom boundary location is defined 
by the continuity relation 
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where K is the thermal diffusivity, k the thermal conductivity, 
p the density of the solid, L the latent heat of fusion, f the 
fraction of material melted, T the solidus temperature, and 

Hjj(t) the heat flow into the base of the lithosphere (a function 
of time) . Since we are not interested in the temperature 
solution for very young ages, horizontal heat flow has been 
ignored. 


A solution to this set of equations was obtained by a 
finite difference method utilizing Lagrangian interpolation 
formulae similar to those described by Crank (1956). The bottjm 
heat flow (H^) largely determines the rate and direction of 
phase boundary movements. Taking * 0 , we could check our 
results with an analytic solution (Carslaw and Jaeger, 1959). 

Depth intervals of 1 km and time intervals of 2500 years vrere 
found to be adequate to reproduce the analytic result. For 
H^>0 and constant in time all solutions were found to asymptotically 
approach a constant surface heat flow given by q s (t*®) = kT a /Z(“) = H^. 


To simulate the heat advected to the base of the lithosphere 
by a plume we assumed a Gaussian shaped heat anomaly with a 
center at time t and a half width of 2.34o where t is the time 
when the heat anomaly is a maximum, and o is a variable parameter 
that determines the width of the anomaly. The position of the 
phase boundary and the surface heat flow were calculated for 
a number of bottom heat functions. The first bottom heat 


H, 


By changing the heat function 


H 0 +H p ex P (“• £^1 2 >. we found that the phase boundary 


function used was H^ft) 
tO Hfctt) «o -p— * • -v 0 

would rapidly thin to a maximum thinning at a time less than 5 
my after the heat flow maximum and then gradually thicken at 
later times. Figure 13 shows the position of the phase boundary 
and the surface heat flow as a function of time for the two 
different bottom heating functions. Values used in these models 
were H Q = 0.65 HFU , =3.2 HFU , tp = 90 my and o = 3 my. This 

thinning mechanism produces only a slight perturbation to the 
normal surface heat flow with a maximum deviation in surface 
heat flow occurring about 25 my after the plume is introduced. 


The thinning data on Figure 12 can be approximated by this 
model if several assumptions are made. First, we have assumed 
that the thermal signature of the heat sources is Gaussian. 
Other functions can be assumed. Secondly, since the anomalous 
heating parameters necessary to account for the thinning depend 
upon the "normal" thermal structure of the lithosphere, some 
thermal model for the oceanic lithosphere must be assumed. The 
model assumed is that in Figure 13 with H^lt) = H p = 0.65 IIFU 
and 5% partial melting. The model is somewhat arbitrary, but 
it does match the empirical surface heat flow and topography 
versus age relations (Parsons and Sclater, 1977) . The actual 
lithosphere would be difficult to model in this area because 
the gradient of lithospheric age is almost perpendicular to 
the present direction of plate motion. 


Figure 14 compares the model results with the thinning 
of Figure 12. The rapid thinning is easily accounted for by 
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Figure 14. Sol'd line is the position of the lithosphere/asthenosphere boundary with respect to the 
mean position for the total area. Dotted line is the thinning caused by introduction of 
a heat anomaly at time » 0. Thinning is the same as Fig. 10. Profile C-D was converted 
from distance to time by using a plate velocity of 8.9 cm/yr. H marks the position of 
Hawai i. 
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the model as well as the later thickening. Ages are calculated 
from Hawaii assuming a plate velocity of 8.9 cm/yr (Minster 
et al. , 1974). An interesting but maybe coincidental feature 
of the model is that a good fit is obtained by moving the heat 
flow maximum directly beneath Hawaii. Notice also that the 
maximum thinning occurs about 5 my after the heat flow maximum. 

It is also evident that the simple shape assumed for the heat 
anomaly does not completely reproduce the details of the thinning. 

3 . 4 Concluding Remarks 

We have shown that the portion of the geoid that is co- 
herent and in phase with the topography can be accounted for 
by an Airy isostatic model with a compensation depth of 50-90 
km. The lithospheric thinning predicted by the Airy model is 
consistent with the observed subsidence of the island chain. 

The rapid thinning of the lithosphere under Hawaii predicted by 
the topography is easily accommodated by an increased heat flow 
from the mantle. A thermal model has been developed which 
predicts the movement of a phase boundary due to advected heat. 

Lithospheric thinning and anomalous mantle heat flow may 
account for other areas of anomalous depth than Hawaii. Others 
(Menard, 1973) have proposed that the dynamic forces are res- 
ponsible for some depth anomalies. Local differences in mantle 
heat flow can also explain the same depth anomalies. Haigh 
(1973) suggested that variations in upper mantle temperatures 
might account for depth anomalies in the North Atlantic; this 
is just a slightly different approac’: to our heating anomalies. 

If lithospheric thinning is called upon to explain regions of 
anomalous depth in the oceans, mantle flow patterns need not 
supply large supporting forces. There is currently no evidence 
to suggest that heat flow from the mantle is locally or even 
regionally uniform. 


4. OCEANIC LITHOSPHERIC AND ASTHENOSPHERIC MODELLING 
4 . 1 Analysis 

In recent years, thermomechanical models have been developed 
at UCLA which prescribe the velocity of a rigid lithosphere, and 
analyze the resulting flow and temperature patterns in the 
asthenosphere under different rheological, thermal, and mechanical 
conditions: dependences of viscosity on temperature and stress, 

distribution of heat sources, depth of return flow, etc. 

(Schubert et al., 1976; Froidevaux et al., 1977; Schubert et al., 
1978). The principal observables predicted by these models 
are the topography and the heat flow. In addition, the tempera- 
ture distributions ir. the upper mantle should be reasonably 
consistent with those inferred from seismic Q and V s , and the 
stresses should be reasonably compatible with rock strengths 
measured in the laboratory. 
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In these models, the gravity anomaly Ag is not explicitly 
calculated since for computational convenience the upper boundary 
is taken as rigid. However, the topography is estimated from 
the pressure p' generated by the flow. The pressure p* in turn 
depends on the excess temperature at depth, T-T re f For the 

depth D cf the ocean floor, with respect to a reference level 
D r . (Schubert et al . , 1976, eq. 4; Schubert et al., 1978, 
eqST 27-29) s 

-(p'-p* ref ) 

D_D ref * g (p-p gw ) 


pa 

^ P " P SW ) 



mmCn 


(T“ T ref ) dy 


(4.1) 


where a is the coefficient of thermal expansion. The gravity 
anomaly in turn will result from a combination of the mas3 
excess D-D re f constituted by the topography plus the mass 
deficiency constituted by the density of thermal expansion 
-pa(T-T r f ) at depth. For a wave number k (e.g., Garland, 
1971, p. 173): 


Ag. - 2IIG[-/ e' ky 6p. (y)dy 
K D 

max 


+ / e“ ky ap <$T k (y)dy] (4.2) 

.00 

where Sp^ and 61^ are the Fourier transforms with respect to x (the horizontal 
direction) of P(y) - Psw 31x1 T-T re f respectively. The wave number k corresponds 
to a spherical harmonic degree l of ak, where a is the earth '3 
radius. Hence to obtain the disturbing potential and geoid 

height N^: 

Ag^ = (ak-l)6v./a (4.3) 

and 

N k - aAg k /g(ak-l) 5 Ag k /kg (4.4) 

To integrate eq. 4.2 for a specific flow model requires 

specification of its wave number spectrum k and for each wave 

number the corresponding contribution 6T k to the anomalous 

temperature T-T re f, as a function of depth y. The models in 

Schubert et al. (1976, Fig. 7) for lOcm/yr spreading velocity suggest 

a form for T-T , of: 

ref 

6T * 200sin ^QftO (4.5) 

for t in My and y in km, y<120 km. The wave length \ is 
about 104 km, and hence the wave number k = 2 H/ X is about 
6xl0"4 km“l. The resulting gravity anomaly from eq. 4.2 is 
thus about +13 mgals, and the corresponding geoid height about 
22 meters. 
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A mo> e accurate result would be obtained by numerically 
integrating eq. 4.2, of course. Furthermore, the full solution 
of the problem requires consideration of the effective lower 
boundary, as emphasized by McKenzie (1977) . It is planned to 
explore further the implications of different spreading rates, 
heat source distributions, and rheologies for the gravity field. 

4 . 2 Average Oceanic Lithosphere 

Preliminary results have been obtained for lithospheric 
admittances by integrating the thermal structure of the simplest 
tl/2 lithospheric model (Oldenburg, 1975). The method, however, 
was found to have several shortcomings: 1) Extreme sensitivity 

to the seafloor topography requires that the crustal thermal 
structure be specified at close (50 m) vertical intervals. 2) 
Admittances for wavelengths shorter than 500 km are dependent 
upon the details of the ridge crest model. 3) Stability of the 
admittances is, of course, best for smaller horizontal grid 
spacings. Numerical results for thermal models are commonly 
only calculated at la-ge horizontal spacings to save money; 
the models must 'hen be interpolated, thus introducing in- 
accuracies. 

1/2 

Since both gravity and bathymetry vary as t except at 
the ridge crest and very far from the ridge crest, it is expected 
that the admittance will be nearly constant with wavelength. 

The slight deviations from t^/2 at old ages will affect the 
largest wavelengths. The extent to which the long wavelengths 
differ from a constant will depend upon the length of the pro- 
file; longer wavelengths should have smaller admittances as 
gravity appears to flatten out more than bathymetry after 80 my. 
If admittances are calculated using the enpirical free-air 
gravity and bathymetry vs. age relations of Cochra.i and Talwani 
(1977) the smaller admittances at longer wavelengths cannot be 
discriminated. The slopes of the t 1 /* lines for gravity and 
bathymetry indicated admittances of about 0.01 mgal/m, however, 
Fourier analysis reveals admittances closer to 0.008 mgals/m 
with a few larger than 0.01 mgal/m. The discrepancy can result 
either from the instability of the Fourier transform or from 
the gross averaging of the data, the former being the more 
likely. In either case, the value of 0.007-0.01 mgal/m is 
quite comparable to that obtained from a typical lithosphere 
model (Oldenburg, 1975); however, no wavelength variation can 
be reliably discerned. 

The application of admittances to the thermal evolution of 
the lithosphere will prove to be beneficial only if reliable 
data can be obtained for comparison. Few places exist in the 
oceans where continuous profiles can be obtained from a ridge 
crest to very old lithosphere and it is this very old lithosphere 
that contains the information regarding deviations from the t*/2 
variations. If a profile is to be pieced together, care must 
be taken to assure that no absolute changes in gravity or 
bathymetry are incorporated into the profile. 
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